rm(list = ls())
## llamando a las librerias que (creemos que) vamos a usar
options("install.lock"=FALSE)
if (!require ("pacman")) install.packages ("pacman")

pacman::p_load (dplyr
                , haven
                , ggplot2
                , ggpubr
                , ggthemes  ## si quieren mas themes
                , readstata13
                , readxl
                , sf
                , tidyverse
                , tidyr
                , units
                , viridis ## paleta de colores Viridis
                , wesanderson## p/usar paleta de colores de Wes Anderson
                , stringr
                , RColorBrewer
                , patchwork
                , Rmisc
                , lfe
                , stargazer
                , AER
                , haven
                , skimr
                , modelsummary
                , terra
                , fixest
                , vtable
                , did
                , cowplot
                , grid)

setwd("/Users/mateovasquez/Dropbox/")

#####

base <- read.dta13("master_violence_landtitle_caracteristicas.dta")

#Outcome Variable
names(base)

base <- base %>% 
  dplyr::mutate(total_violence = rowSums(across(c(sumattacks_public, sumattacks_guerrilla,sumattacks_paramilitary)), na.rm = TRUE),
                health_inst = rowSums(across(c(health_post, health_center,hospital)), na.rm = TRUE))

#base <- base %>%
#  mutate(total_violence = rowSums(select(., starts_with("sumattacks_")), na.rm = TRUE),
#         health_inst = rowSums(select(., starts_with("health_")), na.rm = TRUE))

#Renaming Variables 

base <- base %>% 
  dplyr::rename("indian_title_area_cs" = "area_resguardos_indigenas_cs",
                "indian_title_area" = "area_resguardos_indigenas",
                "afro_title_area_cs" = "area_titul_comunidades_negras_cs",
                "afro_title_area" = "area_titul_comunidades_negras")

#Treatment 	

base <- base %>% 
  dplyr::mutate(post = ifelse(Año>1995,1,0),
                any_titleXpost = any_title*post,
                any_title2 = any_title,
                post2 = post)


#the mean level of violence in neighboring municipalities (1988–1995)
Neighbors <- base %>% 
  dplyr::filter(Año<=1995) %>% 
  dplyr::group_by(codmpio) %>% 
  dplyr::summarise(spatlagguerrsum100_100_1995 = mean(spatlagguerrsum100_100, na.rm = T))

summary(Neighbors$spatlagguerrsum100_100_1995)

base <- left_join(base,Neighbors,
                  by = "codmpio")

#Controls X post
base <- base %>%
  dplyr::mutate(
    altura_post = altura*post, #ELEVATION
    rainfall_post = rainfall*post, #RAINFALL
    coca_0_1_post = coca_0_1*post, # COCA
    presence_mines1560_post = presence_mines1560*post, #GOLD MINES
    num_encomiendas1560_post = num_encomiendas1560*post, #ENCOMIENDAS
    royalroaddist_post = royalroaddist*post, #distance to royal roads
    distancia_mercado_post = distancia_mercado*post, #nearest market towns (km)
    pop95_post = pop95*post, #municipal population in 1995 
    pct_minority85_post = pct_minority85*post, #share of minority population in 1985 (%)
    conflictos_1901_1917_post = conflictos_1901_1917*post,#historical conflict using data on the inci- dence of land conflict between 1901 and 1931
    conflictos_1918_1931_post = conflictos_1918_1931*post,#historical conflict using data on the inci- dence of land conflict between 1901 and 1931
    Violencia_48_a_53_post = Violencia_48_a_53*post, #presence of armed combat during La Violencia (1948–1953)
    anucraids1971_78_post = anucraids1971_78*post, #frequency of land invasions between 1971 and 1978
    spatlagguerrsum100_100_1995_post = spatlagguerrsum100_100_1995*post, #the mean level of violence in neighboring municipalities (1988–1995)
    ltotal_plots_rfmd_1960_85_post = ltotal_plots_rfmd_1960_85*post, # total number of plots reformed between 1960 and 1985 to proxy for rural grievances.
    totviolencia85_post = totviolencia85*post) #ojo: La Violencia (1948 - 1958)

# Applying asinh transformation to each specified variable directly in mutate

base <- base %>%
  mutate(
    as_total_violence = asinh(total_violence),
    as_sumattacks_public = asinh(sumattacks_public),
    as_sumattacks_paramilitary = asinh(sumattacks_paramilitary),
    as_sumattacks_guerrilla = asinh(sumattacks_guerrilla))

# Running the regression with felm, including fixed effects and clustering

base <- base %>%
  mutate(Depto_year = paste(coddepto, Año, sep = "_"))

m1 <- felm(as_sumattacks_paramilitary ~ any_titleXpost:tax_office +  
                 any_titleXpost +
                 tax_office +
                 any_title2 + post2 +
             altura_post + 
             altura + 
             rainfall_post + 
             rainfall + 
             coca_0_1_post + 
             coca_0_1 + 
             presence_mines1560_post + 
             presence_mines1560 + 
             num_encomiendas1560_post  + 
             num_encomiendas1560 + 
             royalroaddist_post + 
             royalroaddist + 
             distancia_mercado_post + 
             distancia_mercado + 
             pop95_post + 
             pop95 + 
             pct_minority85_post + 
             pct_minority85 + 
             conflictos_1901_1917_post +
             conflictos_1901_1917 +
             conflictos_1918_1931_post +
             conflictos_1918_1931 +
             anucraids1971_78_post + 
             anucraids1971_78 + 
             totviolencia85_post +
             totviolencia85 +
             spatlagguerrsum100_100_1995_post + 
             spatlagguerrsum100_100_1995 + 
             ltotal_plots_rfmd_1960_85_post +
             ltotal_plots_rfmd_1960_85 | codmpio + Año | 0 | codmpio, data = base)

m2<- felm(as_sumattacks_paramilitary ~ any_titleXpost:pave_priroads95 +  
            any_titleXpost +
            pave_priroads95 +
            any_title2 + post2 +
            altura_post + 
            altura + 
            rainfall_post + 
            rainfall + 
            coca_0_1_post + 
            coca_0_1 + 
            presence_mines1560_post + 
            presence_mines1560 + 
            num_encomiendas1560_post  + 
            num_encomiendas1560 + 
            royalroaddist_post + 
            royalroaddist + 
            distancia_mercado_post + 
            distancia_mercado + 
            pop95_post + 
            pop95 + 
            pct_minority85_post + 
            pct_minority85 + 
            conflictos_1901_1917_post +
            conflictos_1901_1917 +
            conflictos_1918_1931_post +
            conflictos_1918_1931 +
            anucraids1971_78_post + 
            anucraids1971_78 + 
            totviolencia85_post +
            totviolencia85 +
            spatlagguerrsum100_100_1995_post + 
            spatlagguerrsum100_100_1995 + 
            ltotal_plots_rfmd_1960_85_post +
            ltotal_plots_rfmd_1960_85 | codmpio + Año | 0 | codmpio, data = base)

m3<- felm(as_sumattacks_paramilitary ~ any_titleXpost:agbank_office +  
            any_titleXpost +
            agbank_office +
            any_title2 + post2 +
            altura_post + 
            altura + 
            rainfall_post + 
            rainfall + 
            coca_0_1_post + 
            coca_0_1 + 
            presence_mines1560_post + 
            presence_mines1560 + 
            num_encomiendas1560_post  + 
            num_encomiendas1560 + 
            royalroaddist_post + 
            royalroaddist + 
            distancia_mercado_post + 
            distancia_mercado + 
            pop95_post + 
            pop95 + 
            pct_minority85_post + 
            pct_minority85 + 
            conflictos_1901_1917_post +
            conflictos_1901_1917 +
            conflictos_1918_1931_post +
            conflictos_1918_1931 +
            anucraids1971_78_post + 
            anucraids1971_78 + 
            totviolencia85_post +
            totviolencia85 +
            spatlagguerrsum100_100_1995_post + 
            spatlagguerrsum100_100_1995 + 
            ltotal_plots_rfmd_1960_85_post +
            ltotal_plots_rfmd_1960_85 | codmpio + Año | 0 | codmpio, data = base)

m4<- felm(as_sumattacks_paramilitary ~ any_titleXpost:hospital +  
            any_titleXpost +
            hospital +
            any_title2 + post2 +altura_post + 
            altura + 
            rainfall_post + 
            rainfall + 
            coca_0_1_post + 
            coca_0_1 + 
            presence_mines1560_post + 
            presence_mines1560 + 
            num_encomiendas1560_post  + 
            num_encomiendas1560 + 
            royalroaddist_post + 
            royalroaddist + 
            distancia_mercado_post + 
            distancia_mercado + 
            pop95_post + 
            pop95 + 
            pct_minority85_post + 
            pct_minority85 + 
            conflictos_1901_1917_post +
            conflictos_1901_1917 +
            conflictos_1918_1931_post +
            conflictos_1918_1931 +
            anucraids1971_78_post + 
            anucraids1971_78 + 
            totviolencia85_post +
            totviolencia85 +
            spatlagguerrsum100_100_1995_post + 
            spatlagguerrsum100_100_1995 + 
            ltotal_plots_rfmd_1960_85_post +
            ltotal_plots_rfmd_1960_85  | codmpio + Año | 0 | codmpio, data = base)

# Calculate means in advance to avoid complexities within the stargazer call
mean_paramilitary_full <- round(mean(base$sumattacks_paramilitary, na.rm = TRUE), 3)

stargazer(m1, m2, m3, m4,
          type = "text",
          title = "Non-Coercive State Capacity",
          keep = c(":"),
          align = TRUE,
          omit.stat = c("rsq", "ser"),
          dep.var.labels = c("Paramilitary Attacks"),
          add.lines = list(
            c("Outcome Mean", mean_paramilitary_full, 
              mean_paramilitary_full, 
              mean_paramilitary_full,
              mean_paramilitary_full)),
          notes = c("The unit of analysis is a municipality-year. We apply the inverse hyperbolic sine (IHS) transformation to all the outcomes to correct for the skewed distribution in the number of attacks. The table presents the main results using the assignment of a land title under Law 70 to define the treatment indicator (Law 70 Title). Robust standard errors clustered at the municipality level are reported in parentheses. *** is significant at the 1% level, ** is significant at the 5% level and * is significant at the 10% level. A “Pacific” sample limits the analysis to the Pacific region.")
)

